!******************************************************************************************
!*** author: p. areias
!*** generalized linear solution with sparse technology
!*** can be wrapped in a wider and customized interface
!*** it uses topological ordering to apply multiple point
!*** constraints
!*** the solution is obtained by 
!*** the sparse lu
!*** the usage is as follows:
!*** 
!*** level of qa (chkn), n equiv
!*** [0] - first reading
!*** [1] - first compiling
!*** [2] - second reading
!*** [3] - first debugging (with academic tests)
!*** [4] - testing (with academic tests)
!*** [5] - exception handling
!*** [6] - ensuring that all features are included
!*** [7] - testing (with real tests)
!*** [8] - commenting and spellcheck
!*** [9] - last complete check 
!*** [10] - full verification and discussion with peers
!*** [11] - marked with quality assurance qa/(dd:mm:yyyy)/author
!*** current level of qa=[9]
!*** qa/(12:10:2010)/pareias
!
!*** to start a new problem:
!    call startlinear(ntgl,nno) 
!    where ntgl is the maximum number of dof types per node
!    and nno is the number of nodes. For example, when analyzing a 3d solid mechanics
!    problem with 3 displacement components, ntgl will be 3. for a shell problem with 
!    6 degrees-of-freedom per node, ntgl will be 6, etc.
!    note that some constituents (elements and mpcs) may only have some DOF types activated
!    and therefore the maximum number of activated node should be used
!  
!*** to add an additive constituent:
!*** 
!******************************************************************************************
module linearcore
  use basfun
  private
  public::store,linearsolver,startlinear,mpctreat,newndof,matp,startinsertion
  save

!*** store interfaces
!--------------------
  interface store
     module procedure storeelement,storempc
  end interface store

!*** symbolic and start flag
!---------------------------
  logical::symbo=.true.
  logical::start=.false.
!---------------------------

!*** relevant dimensions
!------------------------------------------------------------------------------------------------
  integer::nno=0     ! number of nodes
  integer::ntgl=0    ! maximum number of dof types per node
  integer::ndof=0    ! number of dof corresponding to the original assignment
  integer::newndof=0 ! number of dof corresponding to the new assignment
  integer::ntyp=0    ! number of types, which is inferred from the input data
!------------------------------------------------------------------------------------------------

!*** mapping between degrees-of-freedom and types
!*** types can be related to the original NODE / TYPE of dof assignment
!------------------------------------------------------------------------------------------------
  integer,dimension(:),allocatable::dofoftype,typeofdof ! dof of type, type of dof
!------------------------------------------------------------------------------------------------

!*** elements and additive constituents
!-------------------------------------------------------------------------------------------------------
  integer::nelm=0                                     ! number of elements
  integer::nelmfake=0                                 ! total number of fake elements, nelm+nmpc 
  integer,dimension(:),allocatable::elmp,elmd,elmdtyp ! element pointer,dof,type
  integer,dimension(:),allocatable::newelmp,newelmd   ! element new pointer, new dof 
! (final numbering used and it will be compressed)
  integer,dimension(:),allocatable::elmmp             ! element stiffness pointer
  real(8),dimension(:),allocatable::elmmat,elmvec     ! element matrix and vector
!-------------------------------------------------------------------------------------------------------

!*** mpcs
!------------------------------------------------------------------------------------------------
  integer::nmpc=0                                     ! number of mpcs
  integer,dimension(:),allocatable::slvd,slvdtyp      ! slave dof,type
  integer,dimension(:),allocatable::masp,masd,masdtyp ! master pointer,dof,type
  integer,dimension(:),allocatable::mpcmp             ! mpc matrix pointer
  real(8),dimension(:),allocatable::mpcmat,mpcvec     ! mpc matrix and vector (both vectors)
  real(8),dimension(:),allocatable::mpcmat2           ! second derivatives - these are matrices
!------------------------------------------------------------------------------------------------

!*** dof destinations and storage
!-------------------------------------------------------------------------------------------------------
  integer,dimension(:),allocatable::destp,destd           ! destination pointer, destination dof
  integer,dimension(:),allocatable::newdestp,newdestd     ! new destination pointer, new destination dof
! (prior to new numbering is used)
  real(8),dimension(:),allocatable::destvec,destmat       ! destination vector, destination matrix
  real(8),dimension(:),allocatable::newdestvec,newdestmat ! new destination vector, new destination matrix
!-------------------------------------------------------------------------------------------------------

!*** dof topological ordering ordering
!------------------------------------------------------------------------------------------------
  integer,dimension(:),allocatable::top
!------------------------------------------------------------------------------------------------

!*** global stiffness
!-------------------------------------------------------------------------------------------------
  integer,dimension(:),allocatable::matp,matd   ! matrix pointer and destination row-by-row
  real(8),dimension(:),allocatable::mat         ! matrix coefficients
  integer,dimension(:),allocatable::ilin,icol   ! row and column format for certain sparse solvers
!-------------------------------------------------------------------------------------------------

!*** reduced right hand side (overlapped with solution)
!------------------------------------------------------------------------------------------------
  real(8),dimension(:),allocatable::vec 
!------------------------------------------------------------------------------------------------

!*** new dof numbers
!------------------------------------------------------------------------------------------------
  integer,dimension(:),allocatable::newdof,olddof
!------------------------------------------------------------------------------------------------

!*** global stiffness addressing by the cliques
!------------------------------------------------------------------------------------------------
  integer,dimension(:),allocatable::cliquep,cliqued
!------------------------------------------------------------------------------------------------

!*** force residual
!--------------------------
  real(8)::residforc=0.0d00
!--------------------------

!*** iempc is the number of mpcs inserted
!----------------------------------------
  integer::iempc=0
!----------------------------------------
contains

!***************************************************
!*** warns that the problem was not properly started
!***************************************************
  subroutine warnnotstart
    write(*,*)"prior to use of linearcore, the subroutine"
    write(*,*)"startlinear must be invoked"
    stop	
  end subroutine warnnotstart
  subroutine type_of_dof(ieq,ino,ity)
    jeq=typeofdof(olddof(ieq))
    ino=1+(jeq-1)/ntgl
    ity=ipermut(ntgl,jeq)
  end subroutine type_of_dof
!***********************
!*** topological sorting
!*** c dieter jungnickel
!*** chk0 
!*** rpt
!***********************
  subroutine doftop(na,ai,aj,acyclic,top)
    integer,dimension(*)::ai,aj
    integer,dimension(:),allocatable::ind,l,top
    logical::acyclic
    m=1
    call salloc(na,ind)
    call salloc(na,l)
    do i=1,na
       if(aj(ai(i)).ne.i)then
          do j=ai(i),ai(i+1)-1
             if(aj(j).gt.0)ind(aj(j))=ind(aj(j))+1
          end do
       end if
    end do
    ik=0
    do i=1,na
       if(ind(i).eq.0)then
          ik=ik+1
          l(na+1-ik)=i
       end if
    end do
    mk=na
    call salloc(na,top)
    do while(ik.ne.0)
       i=l(mk)
       mk=mk-1
       ik=ik-1
       top(m)=i
       m=m+1
       if(aj(ai(i)).ne.i)then
          do j=ai(i),ai(i+1)-1
             ig=aj(j)
             if(ig.gt.0)then
                ind(ig)=ind(ig)-1
                if(ind(ig).eq.0)then
                   ik=ik+1
                   l(mk+1-ik)=ig
                end if
             end if
          end do
       end if
    end do
    if(m.eq.na+1)then
       acyclic=.true.
    else
       acyclic=.false.
    end if
    do i=1,na/2
       i1=top(na+1-i)
       top(na+1-i)=top(i)
       top(i)=i1              
    end do
    deallocate(ind,l)
  end subroutine doftop

!*******************************************
!*** address of a many-to-many relation
!*** j=0 requests the number of destinations
!*** chk0
!*******************************************
  integer function mmaddress(ipo,i,j)
    integer,dimension(*)::ipo
    if(j.eq.0)then
       mmaddress=ipo(i+1)-ipo(i)
    else
       mmaddress=ipo(i)-1+j
    end if
  end function mmaddress

!***********************
!*** clique addressing 
!*** find all the matrix
!*** pointers
!*** chk0
!***********************
  subroutine clqaddress(nel,ind,jnd,istif)
    integer::nel
    integer,dimension(:)::ind,jnd
    integer,dimension(:),allocatable::istif
    call salloc(nel+1,istif)
    do iel=1,nel
       itemp0=ind(iel)
       itemp1=ind(iel+1)
       jtemp0=jnd(iel)
       jtemp1=jnd(iel+1)
       istif(iel)=(itemp1-itemp0)*(jtemp1-jtemp0)
    end do
    call mudlis(nel,istif)
  end subroutine clqaddress

!****************************
!*** clique addressing 
!*** locates the position
!*** of a clique component
!*** in a large matrix
!*** composed of many cliques
!*** chk0
!****************************
  integer function indstiff(istif,ind,iel,iedof,jedof)
    integer,dimension(:)::istif,ind
    indstiff=istif(iel)-1+id2d(mmaddress(ind,iel,0),iedof,jedof)
  end function indstiff

!********************************************************
!*** numeric assembling (updates matrix)
!*** the clique addresses clqp,clqd indicate the position
!*** chk0
!********************************************************
  subroutine nmassb(iel,ind,clqp,clqd,estif,matrix)
    implicit real(8)(a-h,o-z)
    integer,dimension(:)::ind,clqp,clqd
    real(8),dimension(:)::estif,matrix
    nedof=ind(iel+1)-ind(iel)
    do jedof=1,nedof
       do iedof=1,nedof
          iz=clqd(indstiff(clqp,ind,iel,iedof,jedof))
          matrix(iz)=matrix(iz)+estif(id2d(nedof,iedof,jedof))
       end do
    end do
  end subroutine nmassb

!***********************************************
!*** symbolic assembling 
!*** will provide clqp and clqd from lep and led
!*** chk0
!***********************************************
  subroutine symassemb(nel,lep,led,clqp,clqd,neq,mp,md)
    integer,dimension(:)::lep,led
    integer,dimension(:),allocatable::iw,letp,letd,ijle,mp,md,clqd,clqp
    if(neq.le.0)then
       write(*,'(a)') "symbolic assembling only works with a-priori established number of dofs"
       stop
    end if
    if(nel.le.0)then
       write(*,'(a)') "symbolic assembling only works with more than 0 elements"
       stop
    end if
    call invind(nel,neq,lep,letp,led,letd,ijle)
    call clqaddress(nel,lep,lep,clqp)
    call atimesb1(neq,letp,letd,nel,lep,led,igash,mp)
    call salloc(neq,mp,md)
    call salloc(numinj(nel,lep,led),iw)   
    l=0
    do ira=1,neq
       do iza=letp(ira),letp(ira+1)-1
          iel=letd(iza)
          igl=ijle(iza) 
          jgl=0
          do izb=lep(iel),lep(iel+1)-1
             jgl=jgl+1
             mpb=led(izb) 
             ip=iw(mpb)
             if(ip.eq.0)then
                l=l+1
                md(l)=mpb
                iw(mpb)=l
                llp=indstiff(clqp,lep,iel,igl,jgl)
                call insert(clqd,llp,l)
             else
                llp=indstiff(clqp,lep,iel,igl,jgl)
                call insert(clqd,llp,ip)
             end if
          end do
       end do
       do izc=mp(ira),l
          iw(md(izc))=0
       end do
    end do
    call salloc(0,iw)
  end subroutine symassemb

!**************************************************************************************
!*** expands constrained solution into unconstrained
!*** to obtain the global solution
!*** chk0  
!**************************************************************************************
  subroutine expands(nodofn,ntyp,olddf,nwdof,typdf,newdestp,newdestd,newdestvec,newdestmat,vec,soluc)
    use basfun
    implicit real(8) (a-h,o-z)          
    integer,dimension(:)::nwdof,typdf,newdestp,newdestd,olddf
    real(8),dimension(:)::newdestvec,newdestmat,vec,soluc
!*** sets the solution to zero
    call pconsr(ntyp,soluc)
!*** now travels the dof to expand the solution
    do i=1,nodofn ! old dofs
       ityp=typdf(i)
       soluc(ityp)=newdestvec(i)
       do j=newdestp(i),newdestp(i+1)-1
          if(newdestd(j).ne.0)then
             itemp=nwdof(newdestd(j))
             if(itemp.ne.0)then
                soluc(ityp)=soluc(ityp)+newdestmat(j)*vec(itemp)
             end if
          end if
       end do
    end do
  end subroutine expands

!**********************************
!*** sets the symbolic operation on
!*** this will force the next
!*** call to compute the new 
!*** address lists
!**********************************
  subroutine startlinear(ng,nn)
    symbo=.true.
    start=.true.
    ntgl=ng
    nno=nn
    iel=0
    imp=0
  end subroutine startlinear

!********************
!*** starts insertion
!********************  
  subroutine startinsertion
    iempc=0
  end subroutine startinsertion

!**************************************************
!*** insertion of an element
!*** forces are e-i, external forces minus internal
!*** chk0
!**************************************************
  subroutine insertelement(iel,ncoun,ltype,efor,emat)
    implicit real(8) (a-h,o-z)
    integer,dimension(:)::ltype
    real(8),dimension(:)::efor
    real(8),dimension(:,:)::emat
    if(symbo)then
!*** starts elmp if needed
       if(.not.allocated(elmp))then
          allocate(elmp(1))
          elmp(1)=1
       end if
!*** starts elmmp if needed
       if(.not.allocated(elmmp))then
          allocate(elmmp(1))
          elmmp(1)=1
       end if
!*** updates number of elements
       if(iel.eq.1)nelm=0
       nelm=max(nelm,iel)
!*** stores vector information
       istart=elmp(iel)
       call insertcopy(istart,istart-1+ncoun,elmvec,efor)
       call insert(elmp,iel+1,istart+ncoun)
       call insertcopy(istart,istart-1+ncoun,elmdtyp,ltype)
!*** stores matrix information
       iend=elmmp(iel)+(elmp(iel+1)-elmp(iel))*(elmp(iel+1)-elmp(iel))
       call insert(elmmp,iel+1,iend)              
       do jcoun=1,ncoun
          do icoun=1,ncoun
             call insert(elmmat,indstiff(elmmp,elmp,iel,icoun,jcoun),emat(icoun,jcoun))
          end do
       end do
    else
       if(ncoun.ne.0)then
          call copy(ncoun,elmvec(elmp(iel)),efor)
          call copy(ncoun*ncoun,elmmat(indstiff(elmmp,elmp,iel,1,1)),emat)
       endif
    end if
  end subroutine insertelement

!**********************
!*** insertion of a mpc
!*** chk0
!**********************
  subroutine insertmpc(ncoun,mtype,ntype,vect,matrix,matrix2)
    implicit real(8) (a-h,o-z)
    integer,dimension(:)::ntype
    real(8)::vect
    real(8),dimension(:)::matrix
    real(8),dimension(:,:)::matrix2
    iempc=iempc+1
    impc=iempc
    if(symbo)then
!*** checks for head-in-tail situations
       if(ncoun.ne.0)then
          do i=1,ncoun
             if(mtype.eq.ntype(i))then
                stop "cannot insert head-in-tail mpc"
             end if
          end do
       end if
!*** starts masp and mpcmp if needed
       if(.not.allocated(masp))then
          allocate(masp(1))
          masp(1)=1
       end if
       if(.not.allocated(mpcmp))then
          allocate(mpcmp(1))
          mpcmp(1)=1
       end if
!*** updates number of mpcs
       if(impc.eq.1)nmpc=1
       nmpc=max(nmpc,impc)
!*** stores vector information
       call insert(mpcvec,impc,vect)
       call insert(slvdtyp,impc,mtype)
       istart=masp(impc)
       call insert(masp,impc+1,istart+ncoun)
!*** now the matrix stuff
       call insertcopy(istart,istart-1+ncoun,masdtyp,ntype)
       call insertcopy(istart,istart-1+ncoun,mpcmat,matrix)
!*** stores matrix information
       iend=mpcmp(impc)+(masp(impc+1)-masp(impc))*(masp(impc+1)-masp(impc))
       call insert(mpcmp,impc+1,iend)
       do jcoun=1,ncoun
          do icoun=1,ncoun
             call insert(mpcmat2,indstiff(mpcmp,masp,impc,icoun,jcoun),matrix2(icoun,jcoun))
          end do
       end do
    else
       mpcvec(impc)=vect
       if(ncoun.ne.0)call copy(ncoun,mpcmat(masp(impc)),matrix)
       if(ncoun.ne.0)call copy(ncoun*ncoun,mpcmat2(indstiff(mpcmp,masp,impc,1,1)),matrix2)
    end if
  end subroutine insertmpc

!***********************************
!*** wrap to insert
!*** warning: efor sense is switched
!*** and will correspond to external
!*** forces
!***********************************
  subroutine storeelement(iel,ncoun,lnods,ltyps,efor,emat)
    implicit real(8)(a-h,o-z)
    integer,dimension(:)::lnods,ltyps
    integer,dimension(ncoun)::ltype
    real(8),dimension(:)::efor
    real(8),dimension(:,:)::emat
    if(.not.start)then
       call warnnotstart()
    end if
    do i=1,ncoun
       ity=ltyps(i)
       ino=lnods(i)
       ltype(i)=(ino-1)*ntgl+ity
       efor(i)=-efor(i)
    end do
    call insertelement(iel,ncoun,ltype,efor,emat)
  end subroutine storeelement

!*************************************
!*** wrap to insert mpcs
!*** mnods: slave node
!*** mtyps: slave local type
!*** ncoun: number of masters
!*** nnods: list of master nodes
!*** ntyps: list of master local types
!*** rhs
!*************************************
  subroutine storempc(mnods,mtyps,ncoun,nnods,ntyps,rhs,trm,trm2)
    implicit real(8)(a-h,o-z)
    integer::mnods,mtyps
    integer,dimension(ncoun)::nnods,ntyps,ntype
    real(8),dimension(:)::trm
    real(8),dimension(:,:)::trm2
    real(8)::rhs
    if(.not.start)call warnnotstart()
    mtype=(mnods-1)*ntgl+mtyps
    do i=1,ncoun
       ity=ntyps(i)
       ino=nnods(i)
       ntype(i)=(ino-1)*ntgl+ity
    end do
    call insertmpc(ncoun,mtype,ntype,rhs,trm,trm2)
  end subroutine storempc

!******************************************************
!*** mpc treatment
!*** from gr=jacobian of eq with respect to retained
!*** and gs=jacobian of eq with respect to slave
!*** and eq=equation
!*** calculates b and tmpctreat(ndofm,ndofs,deqdm,d
!******************************************************
  subroutine mpctreat(ndofm,ndofs,gr,gs,eq,b,t)
    implicit real(8) (a-h,o-z)
    real(8),dimension(ndofs,*)::gr
    real(8),dimension(ndofs,*)::gs
    real(8),dimension(:,:),allocatable::bi
    real(8),dimension(*)::eq
    real(8),dimension(:),allocatable::b
    real(8),dimension(:,:),allocatable::t
    call salloc(ndofs,b)
    call salloc(ndofs,ndofm,t)
    call salloc(ndofs,ndofs,bi)
    call matinv(ndofs,det,gs,bi,.false.)
    call matvec(ndofs,ndofs,bi,eq,b)
    call escvec(ndofs,b,-1.0d00,b)
    call matmat(ndofs,ndofs,ndofm,bi,gr,t,1)
    call escvec(ndofs*ndofm,t,-1.0d00,t)
    deallocate(bi)
  end subroutine mpctreat

!*****************************************************
!*** calculates the effect of the mpcs in the elements
!*** slave=matrix.master+vector
!*** chk0
!*****************************************************
  subroutine getsridofthempcs()
    implicit real(8) (a-h,o-z)
    logical::acyclic=.false.
    integer::s
    if(symbo)then
!*** determines the number of types 
!*** by going through the constituents
       ntyp=0
       do i=1,elmp(nelm+1)-1
          ntyp=max(ntyp,elmdtyp(i))
       end do
       if(.not.allocated(masp))then
          allocate(masp(1))
          masp(1)=1
       end if
       do i=1,nmpc
          ntyp=max(ntyp,slvdtyp(i))
       enddo
       do i=1,masp(nmpc+1)-1
          ntyp=max(ntyp,masdtyp(i))
       enddo
!*** allocates dof destinations
       call salloc(elmp(nelm+1)-1,elmd)
       call salloc(nmpc,slvd)
       call salloc(masp(nmpc+1)-1,masd)
!*** assigns unique dof numbers
!*** to replace the types employed by
!*** the user
       call salloc(ntyp,dofoftype)
       ndof=0
       do i=1,nelm
          do j=elmp(i),elmp(i+1)-1
             ity=elmdtyp(j)
             ieq=dofoftype(ity)
             if(ieq.eq.0)then
                ndof=ndof+1
                dofoftype(ity)=ndof
                ieq=ndof
             end if
             elmd(j)=ieq
          end do
       end do
       do i=1,nmpc
          ity=slvdtyp(i)
          ieq=dofoftype(ity)
          if(ieq.eq.0)then
             ndof=ndof+1
             dofoftype(ity)=ndof
             ieq=ndof
          end if
          slvd(i)=ieq
       end do
       do i=1,nmpc
          do j=masp(i),masp(i+1)-1 
             ity=masdtyp(j)
             ieq=dofoftype(ity)
             if(ieq.eq.0)then
                ndof=ndof+1
                dofoftype(ity)=ndof
                ieq=ndof
             end if
             masd(j)=ieq
          end do
       end do
!*** dofoftype may be over-dimensioned
       call salloc(ndof,typeofdof)
!*** calculates the inverse relation
       do i=1,ntyp
          if(dofoftype(i).ne.0)typeofdof(dofoftype(i))=i
       end do
!*** converts the mpcs to a manageable format
!*** which is destp/destd
       call salloc(ndof+1,destp)
       call pconsi(ndof,destp,1)
       destp(ndof+1)=0
       do i=1,nmpc
          destp(slvd(i))=masp(i+1)-masp(i)
       end do
       call salloc(ndof,destp,destd)
       do i=1,ndof
          do j=destp(i),destp(i+1)-1
             destd(j)=i
          enddo
       enddo
       do i=1,nmpc
          s=slvd(i)
          im=0
          do k=masp(i),masp(i+1)-1
             im=im+1
             destd(destp(s)-1+im)=masd(k)
          enddo
       enddo
!*** at this point we have destp/destd and the constituents are numbered 
!*** opens the space for the values
       call salloc(ndof,destp,destmat)
       call salloc(ndof,destvec)
!*** performs a topological ordering of the mpcs
       call doftop(ndof,destp,destd,acyclic,top)
       if(.not.acyclic)then
          write(*,*) "the linear problem has no solution since the dof graph contains at least a cycle"
          write(*,*) "check the input data of linearcore"
          stop
       end if
!*** top gives the equation of the given position
       call salloc(ndof+1,newdestp)
       call pconsi(ndof,newdestp,1)
!*** now we will estimate (in excess) the number of destinations 
       do i=1,ndof
          ieq=top(i)
          n=0
          do j=destp(ieq),destp(ieq+1)-1             
             n=n+newdestp(destd(j))
          enddo
          newdestp(ieq)=n
       enddo
!*** now determines the new destination numbers
!*** which may be repeated
       call salloc(ndof,newdestp,newdestd)
!*** opens the space for the new values
       call salloc(ndof,newdestp,newdestmat)
       call salloc(ndof,newdestvec)
!*** when dofs point to themselves, this must be explicitely enforced
       do i=1,ndof
          if(destp(i+1).eq.destp(i)+1)then
             if(destd(destp(i)).eq.i)newdestd(newdestp(i))=destd(destp(i))
          end if
       enddo
!*** ATTENTION: up to here this is correct
       do i=1,ndof
          ieq=top(i)
          n=0
          do j=destp(ieq),destp(ieq+1)-1
             jeq=destd(j)
             do k=newdestp(jeq),newdestp(jeq+1)-1
                n=n+1
                newdestd(newdestp(ieq)-1+n)=newdestd(k)
             enddo
          enddo
       end do
!*** now contracts the problem since dof may have been deleted by the mpc
!*** gives a simple access to renumbering
!*** newdof(old)
!*** solutions are given in terms of the new numbering
!*** so that they must be mapped back
!*** the final number of dofs is called newndof
       call salloc(ndof,newdof)
       call salloc(ndof,olddof)
       newndof=0
       do i=1,ndof 
          do j=newdestp(i),newdestp(i+1)-1
             iold=newdestd(j)
             if(newdof(iold).eq.0)then
                newndof=newndof+1
                newdof(iold)=newndof
                olddof(newndof)=iold
             end if
          end do
       end do
!*** the new number of dofs is now newndof
!*** now we can allocate the global right-hand-side
       call salloc(newndof,vec)
!*** new element pointers are required
!*** since the generalized elements 
!*** have different ranges
!*** fake elements
       nelmfake=nelm+nmpc
       call salloc(nelmfake+1,newelmp)
!*** now determines the newelmp
!*** and the maximum number of element degrees of freedom
!*** note that these are generalized elements
!*** and correspond to contractions and expansions of
!*** the original elements
       do i=1,nelm
          l=0
          do j=elmp(i),elmp(i+1)-1
             iold=elmd(j)
             l=l+newdestp(iold+1)-newdestp(iold)
          end do
          newelmp(i)=l
       end do
       do i=1,nmpc
          l=0
          do j=masp(i),masp(i+1)-1
             iold=masd(j)
             l=l+newdestp(iold+1)-newdestp(iold)
          end do
          newelmp(i+nelm)=l
       end do
       call salloc(nelmfake,newelmp,newelmd)
!*** now the new element destinations (newelmd)
       do i=1,nelm
          l=0
          do j=elmp(i),elmp(i+1)-1
             iold=elmd(j)! e.g. old may depend on 10
             do k=newdestp(iold),newdestp(iold+1)-1
                inew=newdestd(k)
                l=l+1
!*** note that element destinations are corrected by the newdof wrap
                newelmd(newelmp(i)-1+l)=newdof(inew)
             end do
          end do
       end do
       do i=1,nmpc
          l=0
          do j=masp(i),masp(i+1)-1
             iold=masd(j)
             do k=newdestp(iold),newdestp(iold+1)-1
                inew=newdestd(k)
                l=l+1
!*** note that element destinations are corrected by the newdof wrap
                newelmd(newelmp(i+nelm)-1+l)=newdof(inew)
             end do
          end do
       end do
    end if
!*** this is the non-symbolic part
!*** actual numbers are involved
    do i=1,ndof
       destvec(i)=0.0d00
       do j=destp(i),destp(i+1)-1
          destmat(j)=1.0d00
       enddo
    enddo
    do i=1,nmpc
       s=slvd(i)
       im=0
       do k=masp(i),masp(i+1)-1
          im=im+1
          destmat(destp(s)-1+im)=mpcmat(k)
       end do
       destvec(s)=mpcvec(i)
    enddo
    do i=1,ndof
       if(destp(i+1).eq.destp(i)+1)then
          if(destd(destp(i)).eq.i)then
             newdestmat(newdestp(i))=1.0d00
          endif
       end if
    enddo
    do i=1,ndof
       ieq=top(i)
       newdestvec(ieq)=destvec(ieq)
       n=0
       do j=destp(ieq),destp(ieq+1)-1
          jeq=destd(j)
          rconst=destmat(j)
          if(ieq.ne.jeq)newdestvec(ieq)=newdestvec(ieq)+rconst*newdestvec(jeq)
          do k=newdestp(jeq),newdestp(jeq+1)-1
             sconst=rconst*newdestmat(k)
             n=n+1
             newdestmat(newdestp(ieq)+n-1)=sconst
          end do
       end do
    end do
  end subroutine getsridofthempcs

!***********************
!*** wrap for assembling
!***********************
  subroutine linearsolver(solve,soluc_,react_,residforc_)
    implicit real(8)(a-h,o-z)
    logical::solve
    real(8),dimension(:,:),allocatable::soluc_,react_
    real(8),dimension(:),allocatable,save::soluc,react
    if(.not.start)call warnnotstart()
    if(symbo)then
       call salloc(ntgl,nno,soluc_)
       call salloc(ntgl,nno,react_)
    endif
    call assembling(solve,soluc,react)
    residforc_=residforc
    if(solve)then
       do i=1,ntyp
          ino=1+(i-1)/ntgl
          ity=ipermut(ntgl,i)
          soluc_(ity,ino)=soluc(i)
          react_(ity,ino)=-react(i)
       end do
    end if
  end subroutine linearsolver

!*************************
!*** allocates stftemp and
!*************************
  subroutine detstftemp(stftemp)
    implicit real(8)(a-h,o-z)
    real(8),dimension(:),allocatable::stftemp
    mgle=0
    do i=1,nelm
       nedof=elmp(i+1)-elmp(i)
       mgle_local=0
       do i1=1,nedof
          ieq=elmd(elmp(i)-1+i1)
          mgle_local=mgle_local+newdestp(ieq+1)-newdestp(ieq)
       end do
       mgle=max(mgle,mgle_local)
    end do
    do i=1,nmpc
       nedof=masp(i+1)-masp(i)
       mgle_local=0
       do i1=1,nedof
          ieq=masd(masp(i)-1+i1)
          mgle_local=mgle_local+newdestp(ieq+1)-newdestp(ieq)
       end do
       mgle=max(mgle,mgle_local)
    end do
    call salloc(mgle**2,stftemp)
  end subroutine detstftemp

!******************************
!*** assembling with mpcs
!*** but does not solve
!*** will automatically perform
!*** chk0
!******************************
  subroutine assembling(solve,soluc,react)
    implicit real(8)(a-h,o-z)
    logical::solve
    real(8),dimension(:),allocatable::soluc,react
    integer,dimension(:),allocatable::localdof
    real(8),dimension(:),allocatable::stftemp
!*** gets rid of the mpcs 
    call getsridofthempcs()
    if(symbo)then	
!*** performs the symbolic assembling process and allocates
!*** the global stiffness
!*** symassemb also removes duplicates
!*** it alters newelmp and newelmd
       call symassemb(nelmfake,newelmp,newelmd,cliquep,cliqued,newndof,matp,matd)
!*** calculates the transpose of matp, matd
!*** allocates stiffness 
       call salloc(matp(newndof+1)-1,mat)
!*** also determines an alternative representation
!*** required by certain sparse solvers
       call convert(newndof,matp,matd,ilin,icol)
    endif
!*** allocates the solution arrays
    call salloc(ntyp,soluc)
    call salloc(ntyp,react)
!*** sets the global stiffness and the right-hand-side to zero
    call pconsr(matp(newndof+1)-1,mat)
    call pconsr(ntyp,react)
    call pconsr(newndof,vec)
!*** finds maximum range of all elements
    call detstftemp(stftemp)
    call salloc(newndof,localdof)  
!*** main element cycle
    do i=1,nelm
       nedof=elmp(i+1)-elmp(i)
       nedof2=newelmp(i+1)-newelmp(i)
!*** sets element stiffness to zero
       stftemp=0.0d00
       do j=1,nedof2
          localdof(newelmd(newelmp(i)-1+j))=j
       end do
!*** travels the 
!*** element and performs
!*** the required operations
       do i1=1,nedof
          itemp=elmp(i)-1+i1
          ieq=elmd(itemp)
          ity=typeofdof(ieq)
          s=elmvec(itemp)
          react(ity)=react(ity)+s
          do j1=newdestp(ieq),newdestp(ieq+1)-1			
             ieq1_old=newdestd(j1)
             ieq1=newdof(ieq1_old)
             t1=newdestmat(j1)
             vec(ieq1)=vec(ieq1)+t1*s
             do i2=1,nedof
                jtemp=elmp(i)-1+i2
                jeq=elmd(jtemp)
                r=elmmat(indstiff(elmmp,elmp,i,i1,i2))
                do j2=newdestp(jeq),newdestp(jeq+1)-1
                   ieq2_old=newdestd(j2)
                   ieq2=newdof(ieq2_old)
                   t2=newdestmat(j2)
                   ist=id2d(nedof2,localdof(ieq1),localdof(ieq2))
                   stftemp(ist)=t1*t2*r+stftemp(ist)
                end do
                vec(ieq1)=vec(ieq1)-t1*r*newdestvec(jeq) 
             end do
          end do
       end do
       do j=1,nedof2
          localdof(newelmd(newelmp(i)-1+j))=0
       end do
!*** effectively performs the numerical assembling
       call nmassb(i,newelmp,cliquep,cliqued,stftemp,mat)
!*** sets the localdof back to zero 
    end do
!*** main mpc cycle
    do i=1,nmpc
       ity=slvdtyp(i)
       nedof=masp(i+1)-masp(i)
       nedof2=newelmp(i+nelm+1)-newelmp(i+nelm)
!*** sets element stiffness to zero
       stftemp=0.0d00
       do j=1,nedof2
          localdof(newelmd(newelmp(i+nelm)-1+j))=j
       end do
       rforc=-react(ity)
       iloc1=0
       do i1=1,nedof
          itemp=masp(i)-1+i1
          ieq=masd(itemp)
          do j1=newdestp(ieq),newdestp(ieq+1)-1			
             iloc1=iloc1+1
             ieq1_old=newdestd(j1)
             ieq1=newdof(ieq1_old)
             iloc2=0
             do i2=1,nedof
                jtemp=masp(i)-1+i2
                jeq=masd(jtemp)
                r=rforc*mpcmat2(indstiff(mpcmp,masp,i,i1,i2))
                do j2=newdestp(jeq),newdestp(jeq+1)-1
                   iloc2=iloc2+1
                   ieq2_old=newdestd(j2)
                   ieq2=newdof(ieq2_old)
                   t2=newdestmat(j2)
                   ist=id2d(nedof2,localdof(ieq1),localdof(ieq2))
                   stftemp(ist)=t1*t2*r+stftemp(ist)
                end do
             end do
          end do
       end do
       do j=1,nedof2
          localdof(newelmd(newelmp(i+nelm)-1+j))=0
       end do
!*** effectively performs the numerical assembling
       call nmassb(i+nelm,newelmp,cliquep,cliqued,stftemp,mat)
!*** sets the localdof back to zero 
    end do
!*** determines the reactions
    do i=1,ndof
       ieq=ndof-top(i)+1
       ity=typeofdof(ieq)
       do j=destp(ieq),destp(ieq+1)-1             
          jeq=destd(j)
          jty=typeofdof(jeq)
          react(jty)=react(jty)+react(ity)*destmat(j)
       enddo
    end do
!*** now solves the problem using 
!*** a sparse solver of the choice (given as an external subroutine)
    residforc=rnormi(newndof,vec)
    if(solve)then
       call sparsesolver2
!*** expands the solution
       call expands(ndof,ntyp,olddof,newdof,typeofdof,newdestp,newdestd,newdestvec,newdestmat,vec,soluc)
    end if
!*** sets symbo to false
    symbo=.false.
    deallocate(localdof)
  end subroutine assembling
!
!*** dumps force
!
  subroutine dumpsforce(ntgl,ndof,olddof,typedf,vec,itp)
    integer,dimension(*)::typedf,olddof
    real(8),dimension(*)::vec
    iu=iunit()
    if(itp.eq.1)then
       open(iu,file="force.txt",status="unknown")
    else
       open(iu,file="displacement.txt",status="unknown")
    end if
    do i=1,ndof
       ieq=typedf(olddof(i))
       ino=1+(ieq-1)/ntgl
       ity=ipermut(ntgl,ieq)
       write(iu,"(A,2I3,E11.3)")"ino,ity",ino,ity,vec(i)
    end do
    close(iu)
  end subroutine dumpsforce
!
!*** dumps stiffness
!
  subroutine dumpstiffness(ntgl,ndof,nnz,olddof,typedf,ilin,icol,mat)
    integer,dimension(*)::typedf,ilin,icol,olddof
    real(8),dimension(*)::mat
    iu=iunit()
    open(iu,file="stiffness.txt",status="unknown")
    do iz=1,nnz
       ieq=typedf(olddof(ilin(iz)))
       jeq=typedf(olddof(icol(iz)))
       ino=1+(ieq-1)/ntgl
       jno=1+(jeq-1)/ntgl
       ity=ipermut(ntgl,ieq)
       jty=ipermut(ntgl,jeq)
       write(iu,"(A,4I3,E11.3)")"ino,ity,jno,jty",ino,ity,jno,jty,mat(iz)
    end do
    close(iu)
  end subroutine dumpstiffness

!******************************************
!*** this stage corresponds to the solution
!*** of the linear system
!******************************************
  subroutine sparsesolver
    implicit real(8)(a-h,o-z)
    real(8),dimension(:),allocatable::ress
    real(8),dimension(:),allocatable::ind,jnd
    real(8),dimension(:),allocatable::vnd
    call salloc(newndof,ress)
    call copy(newndof,ress,vec)
    call solvesky(newndof,matp,matd,mat,ress,vec)
!!    call atimesb(na,ia,ja,nb,ib,jb,nc,ic,jc,a,b,c)
!!    call solvesystem(.true.,symbo,newndof,matp,matd,mat,vec)
  end subroutine sparsesolver
  subroutine sparsesolver2
    implicit real(8)(a-h,o-z)
    common/ma37dd/u,lp,mp,ldiag
    integer,save::liw,nsteps,mtype,la,iflag
    integer,dimension(:),allocatable,save::iw,ikeep,iw1
    real(8),dimension(:),allocatable::w,w2
    real(8),dimension(:),allocatable::at
    u=1.0d00
    nnz=matp(newndof+1)-1
    if(symbo)then
       liw=5*(2*nnz+4*newndof+1)
       call salloc(liw,iw)
       call salloc(3*newndof,ikeep)
       call salloc(2*newndof,iw1)
       iflag=0
       nsteps=0
       mtype=1
       la=1.4*(60*nnz+40*newndof)
       u=1.0d00
       call ma37ad(newndof,nnz,ilin,icol,iw,liw,ikeep,iw1,nsteps,iflag)
    endif
    call salloc(la,at)
    call copy(nnz,at,mat)
    u=-1.0d00
    call ma37bd(newndof,nnz,ilin,icol,at,la,iw,liw,ikeep,nsteps,maxfrt,iflag)
    allocate(w(maxfrt),w2(newndof))
    u=-1.0d00
    call ma37cd(newndof,at,la,iw,liw,w,maxfrt,vec,iw1,nsteps,w2,mtype)
    deallocate(at,w,w2)
  end subroutine sparsesolver2
  subroutine linslv
    integer,dimension(:),allocatable::parent,lnz,lp,li
    real(8),dimension(:),allocatable::lx,d,ux
    call salloc(newndof,parent)
    call salloc(newndof,lnz)
    call salloc(newndof+1,lp)
    call ldmsymbolic(newndof,matp,matd,lp,parent,lnz)
    nfill=lp(newndof+1)
    if(symbo)write(*,*)"nfill,nnz",nfill,matp(newndof+1)-1
    call salloc(nfill,li)
    call salloc(nfill,lx)
    call salloc(nfill,ux)
    call salloc(newndof,d)
    call ldmnumeric(newndof,matp,matd,mat,lp,parent,lnz,li,lx,ux,d)
    call ldllsolve(newndof,vec,lp,li,ux)
    call ldldsolve(newndof,vec,d)
    call ldlltsolve(newndof,vec,lp,li,lx)
  end subroutine linslv
!
!*** triangular decomposition
!*** father of column j is the smallest i such as l_(ij).ne.0 (i>j)
!*** flag marks the column that visited a certain row
!*** it outputs parent and lnz. lp is an additional quantity obtainable from lnz
!
  subroutine ldmsymbolic(n,ap,ai,lp,parent,lnz)
    integer,dimension(*)::ap,ai,parent,lnz,lp
    integer,dimension(:),allocatable::flag
!*** allocates flag
    call salloc(n,flag) 
!*** go through every column
    do k=1,n
!*** father of k still not known
       parent(k)=0
!*** marks k as visited by k itself
       flag(k)=k
!*** resets the number of non-zeros of column k
       lnz(k)=0
!*** goes through every row
       do ip=ap(k),ap(k+1)-1
!*** i is the row number
          i=ai(ip)
!*** select rows with lower numbers only (l)
          if(i.lt.k)then
             do 
!*** if the row was already visited by the column, exits
                if(flag(i).eq.k)exit
!*** if there is no parent because there are no rows
!*** less than i visited during the column i cycle 
!*** then the parent is obviously k
                if(parent(i).eq.0)parent(i)=k
!*** increases the column counting
                lnz(i)=lnz(i)+1
!*** flags i as visited by k
                flag(i)=k
!*** goes to parent of i, which can be k and therefore it cycles 
                i=parent(i)
             end do
          end if
       end do
    end do
!*** now calculates the pointer
    lp(1)=1
    do k=1,n
       lp(k+1)=lp(k)+lnz(k)
    end do
    call salloc(0,flag)
  end subroutine ldmsymbolic
!
!*** ldlnumeric
!
  subroutine ldmnumeric(n,ap,ai,ax,lp,parent,lnz,li,lx,ux,d)
    use basfun
    implicit real(8)(a-h,o-z)
    integer,dimension(*)::ap,ai
    real(8),dimension(*)::ax
    integer,dimension(*)::lp
    integer,dimension(*)::parent
    integer,dimension(*)::lnz
    integer,dimension(*)::li
    real(8),dimension(*)::lx,d,ux
    real(8),dimension(:),allocatable::y,z
    integer,dimension(:),allocatable::pattern
    integer,dimension(:),allocatable::flag
    real(8)::l_ki
    integer::top,p2
    allocate(y(n),z(n))
    call salloc(n,pattern);pattern=1
    call salloc(n,flag)
    flag=0
    do k=1,n!ok
       y(k)=0.0e00
       z(k)=0.0d00
       top=n+1!ok
       flag(k)=k!ok
       lnz(k)=0!ok
       p2=ap(k+1)
       do ip=ap(k),ap(k+1)-1
          i=ai(ip)
          if(i.le.k)then
             y(i)=y(i)+ax(ip)
             len=1
             do 
                if(flag(i).eq.k)exit
                pattern(len)=i
                len=len+1
                flag(i)=k
                i=parent(i)
             end do
             do while(len.gt.1)
                len=len-1
                top=top-1
                pattern(top)=pattern(len)
             end do
          endif
       end do
       d(k)=y(k)
       y(k)=0.0e00
       do 
          if(top.ge.n+1)exit
          i=pattern(top)
          yi=y(i)
          y(i)=0.0e00
          p2=lp(i)+lnz(i)
          do jp=lp(i),p2-1
             y(li(jp))=y(li(jp))-lx(jp)*yi
          end do
          l_ki=yi/d(i)
          d(k)=d(k)-l_ki*yi
          li(p2)=k
          lx(p2)=l_ki
          lnz(i)=lnz(i)+1
          top=top+1
       end do
    end do
    deallocate(y,z)
  end subroutine ldmnumeric
end module linearcore
